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Inviscid instability of an incompressible flow 
between rotating porous cylinders 
to three-dimensional perturbations 


Konstantin IliE0 and Andrey Morguli^ 


Abstract 

We study the stability of two-dimensional inviscid flows in an annulus between two porous 
cylinders with respect to three-dimensional perturbations. The basic flow is irrotational, and 
both radial and azimuthal components of the velocity are non-zero. The direction of the radial 
flow can be from the inner cylinder to the outer one (the diverging flow) or from the outer 
cylinder to the inner one (the converging flow). It had been shown earlier in Ref. [I] that, 
independent of the direction of the radial flow, the basic flow can be unstable to small two- 
dimensional perturbations. In the present paper, we prove first that purely radial flow is stable 
and that flows with both radial and azimuthal components are always stable to axisymmetric 
perturbations. Then we show that both the diverging and converging flows are unstable with 
respect to non-axisymmetric three-dimensional perturbations provided that the ratio of the az¬ 
imuthal component of the velocity to the radial one is sufficiently large. Neutral curves in the 
space of parameters of the problem are computed and it is demonstrated that for any ratio 
of the radii of the cylinders, the most unstable modes (corresponding to the smallest ratio of 
the azimuthal velocity to the radial one) are the two-dimensional ones. We also consider the 
corresponding viscous stability problem and construct an asymptotic expansion of its solutions 
for large radial Reynolds numbers. We compute the first-order viscous correction to inviscid 
eigenvalues and show that the asymptotic results give a good approximation to the viscous eigen¬ 
values even for moderate values of radial Reynolds number, which indicates that the instability 
may be observed in real flows. 


1 Introduction 

In this paper we continue our study of the instability of a steady inviscid flow in an annulus between 
two permeable rotating circular cylinders that had been started in Ref. [Tj. The basic flow has both 
radial and azimuthal components which are independent of the azimuthal angle 6 and inversely 
proportional to the radial coordinate r of the polar coordinates system (r, 6 ) with the origin at the 
common axis of the cylinders. The direction of the radial flow can be from the inner cylinder to the 
outer one (the diverging flow) or from the outer cylinder to the inner one (the converging flow). It 
had been shown in Ref. [T] that this flow can be unstable to small two-dimensional perturbations 
(in the framework of the inviscid theory) provided that the ratio of the azimuthal component of 
the velocity to its radial component is larger than a certain critical value. This two-dimensional 
instability is oscillatory, and the neutral modes represent azimuthal travelling waves. The main 
aim of the present study is to understand what happens if three-dimensional perturbations are 
allowed. In particular, we are interested to investigate whether the most unstable mode (i.e. the 
mode that becomes unstable first when the azimuthal component of the velocity is increased from 
0) is two-dimensional or not, and to determine the axial wave number of the most unstable mode. 
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The stability of viscous flows between permeable rotating cylinders with a radial flow had been 
studied by many authors [21 Elia El El 1318]. One of the main aims of these papers was to determine 
the effect of the radial flow on the stability of the circular Couette-Taylor flow to axisymmetric 
perturbations, and the general conclusion was that the radial flow affects the stability of the basic 
flow: both a converging radial flow and a sufficiently strong diverging flow have a stabilizing effect on 
the Taylor instability, but when a divergent flow is weak, it has a destabilizing effect nisi- However, 
it was not clear whether a radial flow itself can induce instability for flows which are stable without 
it. This question had been answered affirmatively by Fujita, Morimoto & Okamoto [9| and later by 
Gallet, Doering and Spiegel m who had demonstrated that a particular classes of viscous flows 
between porous rotating cylinders can be unstable to small two-dimensional perturbations. 

Later it had been shown [T] that both converging and diverging flows can be linearly unstable 
in the framework of the inviscid theory and that the instability persists if small viscosity is taken 
into account. In Ref. m, a two-dimensional viscous stability problem had been considered, 
and it had been shown that not only the particular classes of viscous steady flows considered in 
[M E] can be unstable to two-dimensional perturbations, but this is also true for a wide class 
of steady rotationally-symmetric viscous flows (without any restriction on angular velocities of 
the cylinders and for both converging and diverging flows). A further development of the two- 
dimensional theory can be found in a recent paper by Kerswell m where, among other things, the 
effects of compressibility and nonlinearity have been considered. 

In the inviscid theory, a purely azimuthal flow with the velocity inversely proportional to r 
is stable not only to two-dimensional perturbations (see, e.g, [l3]) but also to three-dimensional 
perturbations (this can be deduced from the sufficient condition for stability given by Howard 
&: Gupta [H]). In [T|, it had shown that this stable flow becomes unstable to two-dimensional 
perturbations if a radial flow is added and that the instability occurs for an arbitrarily weak radial 
flow (here ‘weak’ means ‘weak relative to the azimuthal flow’). Moreover, in the limit when the 
ratio of the radial component of the velocity to its azimuthal component is small, the instability 
becomes independent of the only geometric parameter of the problem - the ratio of the radii of 
the cylinders. It had also been observed that if the azimuthal component of the basic flow is zero, 
i.e. the flow is purely radial, then it is stable (to two-dimensional perturbations). These facts 
indicate that the instability mechanism cannot be explained in terms of known instabilities (e.g., 
such as shear flow instability or centrifugal instability). The asymptotic behaviour of the unstable 
eigenmodes for weak radial flow shows that the limit when the radial component of the velocity 
goes to zero is a singular limit of the linear stability problem [Ij. Adding a weak radial flow to 
a purely azimuthal one results in formation of an inviscid boundary layer near the inflow part of 
the boundary, and the new unstable eigenmodes (that were absent in the purely azimuthal flow) 
appear within this boundary layer. 

In the present paper, we examine the effect of three-dimension perturbations on the inviscid 
stability properties of the basic flow. In particular, we rigorously prove that (i) the purely radial 
flow is stable to small three-dimensional perturbations and (ii) the basic flow, in which both the 
radial and azimuthal components of the velocity are nonzero, is always stable to axisymmetric 
perturbations. We also compute neutral curves on the plane of parameters of the problem, which 
demonstrate that, the most unstable mode is always two-dimensional. 

The outline of the paper is as follows. In Section 2, we introduce basic equations and formulate 
the linear stability problem. Section 3 contains a linear inviscid stability analysis of both the 
diverging and converging flows. In Section 4, the effect of viscosity is considered. Discussion of the 
results is presented in Section 5. 
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2 Formulation of the problem 


2.1 Exact equations and basic steady flow 

We consider three-dimensional inviscid incompressible flows in the gap between two concentric 
circular cylinders with radii ri and r 2 (^2 > ri). The cylinders are permeable for the fluid and 
there is a constant volume flux 2ttQ (per unit length along the common axis of the cylinders) of 
the fluid through the gap (the fluid is pumped into the gap at the inner cylinder and taken out at 
the outer one or vice versa). Q will be positive if the direction of the flow is from the inner cylinder 
to the outer one and negative if the flow direction is reversed. Flows with positive and negative 
Q will be referred to as diverging and converging flows respectively. Suppose that ri is taken as a 
length scale, rf/\Q\ as a time scale, \Q\/ri as a scale for the velocity and pQ^lr\ for the pressure 
where p is the fluid density. Then the Euler equations, written in non-dimensional variables, have 
the form 


2 

V V 

Ut + UUr H- Uq + WUz -= —Pr, 

V uv 1 

Vt + UVr + -VB + WUz H -=- P9, 


V 

Wt + UWr + -W0 + WWz = -Pz, 
r 


- {ru)j. + -ve + Wz = 0. 


( 1 ) 

( 2 ) 

(3) 

(4) 


Here {r,9,z) are the polar cylindrical coordinates, u, v and w are the radial, azimuthal and axial 
components of the velocity and p is the pressure. 

It is well-known that if the flow domain is bounded by impermeable walls, then one needs to 
impose the standard boundary condition of no normal velocity at the walls. This will guarantee 
that the the resulting initial boundary value problem for the Euler equations is mathematically 
well-posed. What is less known is that if there is a non-zero flow of the fluid through the boundary, 
then not only the normal velocity must be given at the boundary, but some additional boundary 
conditions at the inflow part of the boundary must also be imposed. What conditions should 
be added is a subtle question and there are several answers that lead to mathematically correct 
initial boundary value problems (see, e.g., [mtts]). We will use the boundary condition for the 
tangent component of the velocity, which at hrst approximation corresponds to the condition at 
a porous cylinder (see m) and for which the corresponding mathematical problem is well-posed 
(e.g., [E]). Another important reason for using this condition is that it is consistent with the 
vanishing viscosity limit for the Navier-Stokes equations (see, e.g., [lailg]). More precisely, the 
solution of the Euler equations with a normal velocity prescribed on the entire boundary of the 
flow domain and with a tangent velocity prescribed on the inflow part of the boundary represent 
the leading order term of the asymptotic expansion of the solution of the corresponding viscous 
problem in which all components of the velocity are given on the entire boundary. 

So, our boundary conditions are 


u\ 1 = /3, u\ = B/a, 

\r=l ’ \r=a ’ 

where a = r 2 /ri and jB = Q/\Q\, and either 


u 1=7, lu 1= 0 

\r=l ' ’ \r=l 

for the diverging flow (/3 = 1) or 

V I = y/a, w\ = 

\ r=a ' / ’ I r=a 


(5) 

( 6 ) 

( 7 ) 
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for the converging flow (/3 = — 1). In Eqs. ([ 6 ]) and ([7]), 7 is the ratio of the azimuthal velocity to 
the radial velocity at the inner cylinder, i.e. 7 = Q.ir\/\Q\ (where fli is angular velocity of the 
inner cylinder). 

The problem, given by Eqs. ©-dSl) and ([ 6 ]) or 0 , has the following simple rotationally- 
symmetric solution: 

u{r,6) = j3/r, v{r,6)=^/r, w = 0. ( 8 ) 

In what follows we will investigate the stability of this steady flow. 


2.2 Linear stability problem 


Consider a small perturbation (u,v,w,p) of the basic flow ([ 8 ]). It is convenient to write the linearised 
equations in the terms of perturbation vorticity 


cj = cui + a;2 eg + ws 


(9) 


where e,., bq and are unit vectors in the radial, azimuthal and axial directions, respectively, and 
where 


1 

uji = - we- Vz, 
r 

002 = Uz - Wr, 

W3 = - {irv)r - ue) 


The linearised equation can be written as 


dt -\—2 ^^—2 1 ~ 

rp pZ. j 

ot H- 7 ^ O 0 -\ - Or -^ I W 2 —- 7 , 


^ ~ ] OJ 3 = 0 . 


( 10 ) 

( 11 ) 

( 12 ) 

(13) 

(14) 

(15) 


We seek a solution of Eqs. (I10l) - (ll5p in the form of the normal mode 

{u, V, w, u;} = Re {n(r), {)(r), tc(r), 

where n G Z and A: G M are the azimuthal and axial wave numbers respectively. On substituting 
this into Eqs. (I13I1 - (I15I) . we can rewrite them as 


where 


h{r) + ^ dr^ (rcji) = 0 , 

(16) 

h(r) + ^dr'^ (y) 

(17) 

h(r) + — W 3 = 0 

r J 

(18) 

h(r)=.+ !^, 

(19) 

IjTX 

chi = — w — iki), 
r 

( 20 ) 

Cj2 = iku — Wr, 

( 21 ) 

1 / .X in ^ 

UJ 3 = - {rv}^ - u. 

( 22 ) 
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Equations (fT 6 ]l - (fT 8 i) should be solved subject to the boundary conditions: 



u\ , = 0, 

lr=l ’ 

(23) 


u\ =0. 

1 r=a 

(24) 

and either 

D1 1 = 0 , rt) 1 1=0 

lr=l ’ lr=l 

(25) 

for the diverging flow or 

hi =0, If) 1 =0 

\r=a ’ lr=a 

(26) 


for the converging flow. Eqnations (|16p ~ (j24p together with either (12511 or (j26ll represent an eigen¬ 
value problem for a. If there is an eigenvalue a such that Re{a) > 0, then the basic flow is 
unstable. If there are no eigenvalues with positive real part and if there are no pertnrbations with 
non-exponential growth (examples of non-exponential growth can be fonnd, e.g., in [ 20 ]), then the 
flow is linearly stable. In the next section we analyse this eigenvalne problem. 


3 Analysis of the eigenvalue problem 

The eigenvalue problem formulated above can be reduced to a problem of finding zeros of a certain 
entire fnnction. We will show this first for the divergent flow. 


3.1 Diverging flow {(3 = 1) 

3.1.1 Dispersion relation 

Bonndary conditions ([2^ and Eq. ([20P imply that 

^iL=i=0. (27) 


Now let 


g{r) ~ ^ 3- in'y Inr, 


(28) 


so that h{r), given by (1191) . can be written as h{r) = g'{r)jr. Then the general solntion of Eq. (|16p 
is 


rui = 


where C is an arbitrary constant. This and Eq. (I27p imply that (7 = 0 and, therefore, wi(r) = 0, 
so that we have the relation 

271 

— w — ikv = 0. (29) 

r 

Now we assume that n 7 ^ 0. The case of n = 0 will be treated separately. Using (I29p to eliminate 
w from the incompressibility condition 


we obtain 


Integration of Eq. (flSP yields 


- {ru)r H- V + ikw = 0 , 




^3 = Cl 


(30) 


(31) 

(32) 
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for an arbitrary constant Ci. Equations (j32[) and (12^ have a consequence that 

inru = r{rv)r — (33) 

Finally, we use (I33p to eliminate u from Eq. m- As a result, we get the equation 


Grr + -Gr- (k^ + %]G = Ci F{r) 

fp V I 


where 

and 


G{r) = rv{r) 


F{r) = - dr (j'^e . 


(34) 

(35) 

(36) 


Equation ([33]) allows us to rewrite boundary conditions 

G(l) =0, 

G'(l) = Gie-5(i)^ 

G'(a) =Giae-5(“). 


(for u and v) in terms of G: 


(37) 

(38) 

(39) 


Equation (IM|l together with boundary conditions ([37|) ~ (f39|) represent an eigenvalue problem for a 
(that enters the problem via g{r)). 

The general solution of Eq. (I34p can be written as 

r 

G(r) = Gi j F{s)[In{kr)Kn{ks) - Iniks)Knikr)]sds + G2lnikr) + GsKnikr). (40) 
1 

Here In{z) and Kn{z) are the modified Bessel functions of the first and second kind; G 2 and G 3 
are arbitrary constants (recall that Gi is also arbitrary). Substitution of the general solution into 
boundary conditions ([37|l and (l38]l results in the following two equations: 

G2ln{k)+G3Kn{k)=0, 

G2 kl'r,{k) + G3 kK{k) = Gie-®(^). 


Solving these for Gi and G 2 , we obtain 

G 2 = Gi Kn{k)e-^^^\ G 3 = -Gi In{k)e-^^^\ 
Here we have used the Wronskian relation (e.g. my- 

Fn{z)Kniz) - Ir^{z)K'^{z) = 


(41) 


(42) 


With the help of (I4ip . we can rewrite Eq. (I40p in the form 

G(r) = Gi I^ F{s) [In{kr)Kn{ks) - In{ks)Kn{kr)] sds + [In{kr)Kn{k) - In{k)Kn{kr)] | . 

Substituting this into boundary condition (1391) . we obtain the dispersion relation 

a 

J F{s)k [F^{ka)Kn{ks) — Iniks)K'^{ka)] sds 
1 

+k [l'^{ka)Knik) - In{k)K'^ika)] = 0. (43) 
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This dispersion relation can be further simplified as follows. Let X be the integral entering the 
dispersion relation. Recalling that T(r) is given by Eq. (|36l) and integrating by parts, we obtain 

a 

X = kj^ds [ll,{ka)Kniks) - In{ks)K'^{ka)] sds 

1 

= [l'^{ka)Kn{ks) - In{ks)K'Jka)] “ 

a 

e-sW [l'^{ka)K'^{ks) - l'^{ks)K'^{ka)\ ds 

1 

= [l'^{ka)Kn{k) - In{k)K{ka)] 

a 

-k"^ J [l'^{ka)K'^{ks) - ll,{ks)K'^ika)] s'^ds. 

1 

Here again we have used the Wronskian relation (1421) . Substitution of the above formula for X into 
(j43p yields the final expression for the dispersion relation: 

a 

D{a,n,k,-f,a) = J %{ks)K’^{ka) - X^{ka)K’^{ks)] ds = 0 . (44) 

1 

It can be shown that in the limit A: —)• 0 this reduces to the dispersion relation of the corresponding 
two-dimensional problem (considered in [T]). 

The dispersion relation ()44l) has been obtained under assumption that n 7 ^ 0. Nevertheless, it 
can be shown that this dispersion relation is also valid for the axisymmetric mode, n = 0 . 

The eigenfunction G{r) associated with the eigenvalue a can be written as 

r 

G{r) = Cik J e-<^*V2-m7ins ^j>^(^ks'^Kn{kr) - In{kr)K’^{ks)] ds, 

1 

while the corresponding formula for H (r) = ru{r) is 

r 

H{r)=C,-r f ^'“^[i;(A:s)iL;(fcr)-i;(fcrX(M] s^ds. 

mj 

1 

3.1.2 General properties of the dispersion relations fl44|l 

It has been mentioned in [22] that certain conclusions about a two-dimensional counterpart of (1440 
can be made using the Polya theorem (see problem 177 of Part V in |23|, see also [24] 1. It turns out 
that this theorem also works for (|44l) . It is shown in Appendix A that, for the purely radial flow 
(7 = 0 ), the dispersion relation (j44|) has no roots with non-negative real part, so that there are no 
growing normal modes for the purely radial diverging flow. The same is true for the axisymmetric 
mode, n = 0 (see Appendix A). So, we can restrict our attention to non-axisymmetric perturbations 
for 7 7 ^ 0 . 

Also, using the fact that I-n{z) = In{z) and K-n{z) = Kn{z) (e.g. m), we deduce from (f44ll 
that 

D{a,n,k,a,'y) = D{a,—n,k,a,'y), (45) 

D{a,n,k,a,'y) = D{a,—n,k,a,—'y) (46) 

where the bar denotes complex conjugation. These relations imply that it suffices to consider only 
positive n and 7 . 
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Critical yfor a= 1.5 



Figure 1: Neutral curves for /? = 1 (diverging flow), a = 1.5 and n = 1,..., 13. The region above 
each curve is where the corresponding mode is unstable. 

3.1.3 Numerical results 

As we already know, for 7 = 0, all eigenvalues lie in the left half-plane of complex variable a. 
Numerical evaluation of (I44p confirms this fact and shows that when 7 increases from 0, some 
eigenvalues move to the right, and there is a critical value 'jcr > 0 of parameter 7 at which one of 
the eigenvalues crosses the imaginary axis, so that 

Re(iT) >0 for 7 > jcr and Re((T) <0 for 7 < 7 ^. 

We have computed neutral curves (Re((T) = 0) on the {k, 7 ) plane for several values of the geometric 
parameter a and for n = 1, 2,... 20. For all a, the neutral curves look qualitatively similar to what 
is shown in Fig. [TJ One can see that the neutral curves for a few modes with low azimuthal wave 
number can be non-monotonic functions of the axial wave number k (e.g., n = 1,2,3 in Fig. [1]). 
However, all other modes are strictly increasing functions of k. Let r(A:) be the critical value of 7 
minimized over n = 1 ,..., 20 : 

r(A;) = min 7 cr(n. A:). 

n 

Functions r(A:) for several values of the geometric parameter a are shown in Fig. [2j This figure 
demonstrates the following three things. First, r(A:), for any value of the geometric parameter a, 
is an increasing function, so that its minimum is attained at A: = 0, i.e. for the two-dimensional 
mode. Thus the mode that becomes unstable first when 7 increases from 0 (we will call it the most 
unstable mode) is two-dimensional. Second, for small to moderate values of A: (A: < 10), function 
r(A:) considerably depends on a: on one hand, it decreases when a is increased and seems to tend 
to a limit for large a; on the other hand, it grows when a tends to 1. Third, r(A:), for any value of a, 
becomes a linear function of k for sufficiently large k. Moreover, this linear asymptote is the same 
for all values of a. It is shown in Appendix B that in the limit of large k and n, more precisely, if 

n = ak and k —>■ 00 , 


where a is a positive constant, then 


min7cr-(a;, k) ~ 2.4671 k. 
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Figure 2 : Critical 7 minimized over azimuthal wave numbers n = 1, 2,... 20 (F = min„ 7 ^) for the 
diverging flow and a = 1.25,1.5,2,4,8. The region above each curve is where the corresponding 
flow is unstable. Circles show the asymptotic stability boundary for large k derived in Appendix 

B. 

This asymptotic result is shown by circles in Fig. [2J Evidently, it is in a good agreement with the 
numerical results. The azimuthal wave number n of the most unstable mode (that, for a fixed k, 
becomes unstable first when 7 is increased from 0) depends on both a and k. The results of the 
numerical calculations of this quantity are shown in Fig. [3j The jumps in n correspond to the 
intersection points of the neutral curves for individual azimuthal modes. Figure [3] indicates that, 
for sufficiently large k, the azimuthal wave number of the most unstable mode, n, is independent 
of a and n ^ k. These facts are employed in Appendix B where the asymptotic behaviour of 
eigenvalues for large k is considered. 

The graphs of functions Re(ff(r)) and Im (H{r)) corresponding to the critical value of 7 for 
a = 2 and n = 4 are shown in Fig. [H Evidently, when the axial wave number increases, the 
eigenfunction becomes more oscillatory and concentrated near the inner cylinder (i.e. at the flow 
inlet). This fact is also used in the investigation of the asymptotic behaviour of eigenvalues for 
large k in Appendix B. 

3.2 Converging flow (/3 = —1) 

3.2.1 Dispersion relation 

An analysis similar to what we did for /3 = 1 results in the following dispersion relation 

a 

Di(cT,n,fc, 7 ,a) = k^ J e-V 2 +m 7 ins _ r^(k)KUks)] ds = 0. (47) 

1 

Again, it can be shown that in the limit A; —)• 0 this reduces to the dispersion relation of the 
corresponding two-dimensional problem (see [I])- 

Similarly to how this was done in Appendix B for the diverging flow, it can be shown that the 
dispersion relation (j47p has no roots with non-negative real part (i) for the purely radial converging 
flow (i.e. for 7 = 0 and for all n) and (ii) for the axisymmetric mode (for n = 0 and for all 7 ). 
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Most unstable mode for a=1.25,1.5, 2, 4, 8 



Figure 3: The azimuthal wave number of the most unstable mode, n, for the diverging flow as a 
function of /c for a = 1.25,1.5, 2,4, 8 . 


(a) k=1 



(b) k=10 


(c) k=20 




Figure 4: Eigenfunction H(r) for neutral modes (7 = jcr) for the diverging flow, a = 2 and n = 4: 
(a) /c = 1, 7 = 5.9483; (b) A: = 10, 7 = 35.7645; (c) A; = 20, 7 = 128.2941. 
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Critical yfor a= 1.5 



Figure 5; Neutral curves for (3 = —1 (converging flow), a = 1.5 and n = 1,... ,12. The region 
above each curve is where the corresponding mode is unstable. 


The dispersion relation (I47|) has the same symmetry properties as its counterpart (j44p for the 
diverging flow: 

T>i((T,n, fc,a, 7 ) = Di{a,-n,k,a,-f), (48) 

Di{cr,n,k,a,'y) = Di{a,-n,k,a,-'y). (49) 

These relations imply that we need to consider only positive n and 7 . 

3.2.2 Numerical results 

Numerical results for the converging flow are similar to those for the diverging flow: for 7 = 0 , 
all eigenvalues lie in the left half-plane of complex variable cr; when 7 increases from 0 , some 
eigenvalues move to the right and cross the imaginary axis. In the case of the converging flow, 
we will use parameter ka instead of k. This is convenient because, to a certain extent, it allows 
us to eliminate the dependence of the results on the geometric parameter a. We have computed 
neutral curves (Re((T) = 0) on the {ka,'y) plane for several values of the geometric parameter a 
and for n = 1, 2,... 20. For all a, the neutral curves look qualitatively similar to what is shown for 
a = 1.5 in Fig. [5l We have found that, at least for a = 1.25,1.5, 2,4 and 8 , the neutral curves for 
all azimuthal modes are increasing functions of k (this differs from the case of the diverging flow 
where neutral curves for some low azimuthal modes can have a local minimum, e.g. for the modes 
with n = 1, 2 in Fig. [T]). Let F(fea) be the critical value of 7 minimized over n = 1 , 2 ,... 20: 

r(A;a) = min7cr(n, A:a). 

n 

Functions r(A:o) for several values of the geometric parameter a are shown in Fig. [H The following 
conclusions can be drawn from this figure. First, r(A:a) is an increasing function for any value of 
the geometric parameter a (at least in the range 1.25 < a < 8 ), so that its minimum is attained 
at /c = 0, i.e. for the two-dimensional mode. So, the mode that becomes unstable first when 7 
increases from 0 is two-dimensional. Second, one can see that, for small to moderate values of ka 
[ka < 10 ), the critical value of 7 depends on a: it decreases when a is increased and seems to tend 
to a limit for large a; and it grows when a decreases. Third, for any values of a, r(A:a) becomes a 
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Figure 6 : Critical 7 minimized over azimuthal wave numbers n = 1,2,... 20 (F = xmnn^cr) for 
the converging flow and a = 1.25,1.5, 2,4, 8 . The region above each curve is where the correspond¬ 
ing flow is unstable. Circles represent the asymptotic stability boundary for large ka derived in 
Appendix B. 


linear function for sufficiently large ka, and this linear asymptote is the same for all values of a. 
We show in Appendix B that if 

n = aka as ka ^ 00 , 

where a is a positive constant, then 

min 7 cr(Q!, fca) ~ 2.4671 ka. 

a 

This asymptotic result is shown by circles in Fig. [ 6 l One can see that it is in a good agreement 
with the numerical results even if ka is not very large. The comparison of Figures [2] and [ 6 ] shows 
that the critical values of 7 for the converging flow is slightly higher than that for the diverging 
flow. In this sense, the converging flow is more stable than the diverging flow. The azimuthal 
wave number n of the most unstable mode (that, for a fixed ka, becomes unstable first when 7 is 
increased from 0) depends on both a and ka. The results of the numerical calculations are shown 
in Fig. [3 The jumps in n correspond to the intersection points of the neutral curves for individual 
azimuthal modes. One can see in Fig. [3 that, for sufficiently large ka, the azimuthal wave number 
of the most unstable mode, n, is independent of o and n ~ ka. These facts are used in Appendix 
B. 

4 Effect of viscosity 

Here our aim is to show that for sufficiently high Reynolds numbers the unstable inviscid modes 
found in the previous section give a good approximation to the corresponding viscous modes. We 
will restrict our analysis to the diverging flow and compute the first viscous correction to the inviscid 
eigenvalues in the limit of high radial Reynolds numbers. 

The most general steady rotationally-symmetric viscous flow between rotating porous cylinders 
is given by (see, e.g., m) 

u = —, V = V{r) = (50) 

r r 
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Most unstable mode for a=1.25,1.5, 2, 4, 8 



Figure 7: The azimuthal wave number of the most unstable mode, n, for the converging flow versus 
ka for a = 1.25,1.5,2,4,8. 


where R = \Q\/i^ is the radial Reynolds number, (with u being the kinematic viscosity), and A and 
B are constants which depend on a, R, 71 = Qirf/Q and 72 = (with fli and Q 2 being the 

angular velocities of the inner and outer cylinders). We will consider only the diverging flow, i.e. 
/? = 1 in Eq. (ISOll . 

It can be shown (see, e.g., [H]) that in the limit of high Reynolds numbers, R^ 1, the azimuthal 
component of the velocity is well approximated by 

E(r) = ^ + ^ + O (R-^) (51) 

r a 

where rj = R(1 — r/a) is the boundary layer variable and function / is defined as 


f{z) = (72 - 7i)e ^ 


(52) 


Note that the second term in (j5ip (the boundary layer term) is non-zero only in the layer of thickness 
0{R~^) near the outer cylinder, so that the flow is well approximated by the first term everywhere 
except for this thin boundary layer. The first term represents the inviscid flow that coincides with 
(jH]) (up to replacing 71 with 7 ). Note also that the single inviscid flow (l 8 |) represents the high- 
Reynolds-number limit of each member of a one-parameter family of viscous flows (parametrised 
by 72 ). 

Consider a small perturbation {u,v,w,p) in the form of the normal mode 


{u,v,w,p} 


Re 


{u(r), {)(r), u;(r),p(r)}e'"‘+*’"^+*^^ 


where u, v, w are perturbations of the radial, azimuthal and axial components of the velocity, p 
is the perturbation pressure, n G Z and /c G M. Substituting this into the linearised Navier-Stokes 
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equations yields the eigenvalue problem: 


cr + 

a + 

cr + 


inV (3 


inV P 


/3 


2V 


H— dr j u -- V = -drP+ — I Lu - - - 


R 


u 2in 


inV P 


P 


in 


+ — dr ]v + ^v + n(r)u =- p + — [ Li) -^ + 


R 


V 2in 


■ u 


1 


H— dr]w = —ikp + — Lw, 
r r / R 


dr (ru) + ini + ikr w = 0, 

11 ( 1 ) = 0 , u{a) = 0 , 'O(l) = 0 , i{a), wl(l) = 0 , w{a) = 0 . 


InEqs. (j53]) ~(j55]l. 


L = 


dr"^ 




V{r) 


(53) 

(54) 

(55) 

(56) 

(57) 


For k = 0, this system reduces to the two-dimensional viscous stability problem whose asymptotic 
behaviour has been studied in ffU- An asymptotic expansion of solutions of Eqs. (|53]l - (l57l) can be 
constructed in almost exactly the same manner and has the form 


a = ao + R ^(Ti + O [R ^) , (58) 

u = nS(r) + R-\u[{r) + ^( 77 )] + O {R-^) (59) 

i = {io(r) -h io{r]) + R~^[i[{r) -|- i\{r])] + O {R~^) (60) 

w = WQ{r) + u)q('>j) + R~^lwi(r) -t- + O {R~^) (61) 

P = fo{r) + pUv) + R~^\Pi{r) + p\{r])] + O {R-'^) (62) 


Here r] = R{1 — r/a) is the boundary layer variable, functions with superscript “r” represent the 
regnlar part of the expansion, and functions with superscript “b” give us boundary layer corrections 
to the regular part. The boundary layer part of the expansion exponentially decays outside the 
boundary layer. A brief account of the details of the asymptotic expansion is given in Appendix 
C. Here we will only present the results. 

In Eq. (f58]l . ao is the inviscid eigenvalue discussed in the previous section, and cJi is the 
first-order viscous correction, computed in Appendix C. The exact eigenvalue problem, given by 
Eqs. (I53l) - (j57p was solved numerically using an adapted version of a Fourier-Chebyshev Petrov- 
Galerkin spectral method described in [25]. We have computed the eigenvalue with largest real 
part, a, numerically for a range of values of the Reynolds number R and compared the resnlts 
with the inviscid eigenvalue uq and the first order viscous approximation ao + aijR. The error of 
approximating a by these is shown in Figs. [ 8 ] and [9] where Eo = \a — ao\ and Ei = \a — ao — ai/R\. 
Figure [ 8 ] shows the plots of Eo and Ei versus ii for a = 1.5 and a = 2 for various values of 72 . 
In both plots, n = 1, A: = 1 and 71 = 20. It is evident that in both cases, ao + ai/R gives much 
better approximation for a if the Reynolds number is sufficiently high. One can also see that Ei 
is quite small even for R = 200, which is not a very high value of the Reynolds number. Figure [9] 
shows Eo and Ei as functions of R for o = 2 , n = 1 , 71 = 20 and 72 = 0 and for various values of 
the axial wave number, k. One can see that the dependence of Eo and Ei on k is very weak. In 
fact, all curves shown in Fig. ([9]) almost the same as the curve corresponding to two-dimensional 
perturbation and shown in Fig. 2(b) of Ref. [11] . 

Both Fig. [8| and [9| show that even for moderate values of R such as R = 200 the asymptotic 
formula a ra ao + ai/R (where ai is computed in Appendix C) produces eigenvalues that are very 
close to the eigenvalues of problem (I53p - (I57|) . This means not only that the inviscid instability 
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Figure 8 : E{) = \a — ao] and Ei = \a — ao — cri/R\ versus R. In both plots, n = 1, = 1 and 

7 i = 20. (a) corresponds to a = 1.5 and (b) shows the results for a = 2. 
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Figure 9: Eq = \a — ao\ and Ei = \a — ao — ai/R\ versus R: a = 2, n = 1, 71 = 20, 72 = 0 and 
k = 0.5,1, 2 and 4. 
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studied here persists if viscosity is taken into account, but also that the asymptotic theory works 
well for Reynolds numbers which are not very high, and this, in turn, implies that the instability 
may be observed in real flows. 

5 Discussion 

We have shown that, in the framework of the inviscid theory, a simple rotationally-symmetric flow 
between two permeable cylinders is unstable to small three-dimensional perturbations. We gave 
a rigorous proof of the facts that the purely radial diverging and converging flows are stable and 
that unstable modes cannot be axisymmetric. Numerical calculations demonstrated that (i) for all 
values of the geometric parameter a in the range from 1.25 to 8 , the most unstable mode (i.e. the 
mode that becomes unstable first when parameter 7 is increased from 0 ) is two-dimensional and 
(ii) the critical value of 7 minimized over all azimuthal modes is a strictly increasing function of 
the axial wave number. 

We have also derived an asymptotic expansion of the corresponding viscous stability problem 
for R 3> 1 , computed the first-order viscous correction to inviscid eigenvalues and compared the 
asymptotic results with numerically obtained viscous eigenvalues. This demonstrated that the 
asymptotic results give a very good approximation even for Reynolds numbers that are not partic¬ 
ularly high, such as R = 200, which suggests that the instability may be be observed in real flows. 
Of course, precise conditions under which this instability will be dominating require a further study, 
and this is a subject of a continuing investigation. 

R is known that a purely azimuthal flow with the velocity inversely proportional to r is stable to 
three-dimensional perturbations (this follows from a sufficient condition for stability given in [14j). 
The present paper shows that a purely radial flow is also stable to three-dimensional perturbations. 
These facts indicate that the physical mechanism of the instability must rely on some destabilising 
effect arising from the presence of both the radial and azimuthal components of the basic flow. 
It has been shown in our previous paper [T] that if a small radial component is added to the 
purely azimuthal flow, it immediately becomes unstable for any value of the ratio of the radii of 
the cylinders, and the growth rate is proportional to the square root of the ratio of the radial 
component of the velocity to the azimuthal one. The asymptotic behaviour of two-dimensional 
unstable eigenmodes in the limit of weak radial flow (see m) shows that this limit is a singular 
limit of the linear stability problem. Adding a weak radial flow to a purely azimuthal one results 
in formation of an inviscid boundary layer near the inflow part of the boundary, and the unstable 
eigenmodes are concentrated within this boundary layer. These facts suggest the following physical 
mechanism of the instability: in a purely azimuthal flow there are no unstable eigenmodes, but 
when we add a weak radial flow, this leads to appearance of new unstable eigenmodes (which do 
not exist at all if there is no radial flow) concentrated within a thin inviscid boundary layer near 
the inflow part of the boundary. This mechanism bears some resemblance to the tearing instability 
in the magnetohydrodynamics, where the unstable tearing mode appears when a small resistivity 
is taken into consideration (see Ref. [26]). 

The instability considered here is oscillatory. The two-dimensional neutral modes represent 
azimuthal travelling waves, while the three-dimensional ones are helical waves. An oscillatory 
instability and appearance of azimuthal and helical waves are also present in the Couette-Taylor 
flow between impermeable cylinders. In the Couette-Taylor flow, these waves are observed at 
moderate azimuthal Reynolds numbers and are associated with viscous effects (see, e.g., [27|). The 
results of the present paper show that, in the presence of a radial flow, azimuthal and helical 
waves may appear at arbitrarily large radial Reynolds numbers, which means that these waves can 
be generated not only by fluid viscosity but also by a radial flow. This has a certain similarity 
with self-oscillations observed in numerical simulations of inviscid flows through a channel of finite 
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length [28]. A more detailed analysis of the effect a radial flow on the stability characteristics 
of the Couette-Taylor flow requires a further investigation which would take full account of the 
viscosity. A particularly interesting question that arises in this context is the relation between the 
instability studied here and the classical centrifugal instability that leads to the formation of the 
Taylor vortices. Here is an interesting paradox; in the inviscid theory, axisymmetric modes cannot 
be unstable, but it is well known that the monotonic instability with respect to axisymmetric 
perturbation occurs in the Couette-Taylor flow with radial flow (see, e.g., m)- Our hypothesis is 
that the monotonic axisymmetric and oscillatory non-axisymmetric instabilities are well separated 
in the space of parameters of the problem. If this were so, it would mean that our instability can 
be observed experimentally. This, however, requires a further theoretical study and is a topic of a 
continuing investigation. 

The results presented here are mainly of theoretical interest. However, as was argued in Ref. 
uni, they may be relevant to astrophysical flows such as accretion discs (see also Refs. mm)- 
Our results may also shed some light on the physical mechanism of the formation of strong rotating 
jets in flows produced by a rotating disk which had been observed experimentally (see [SolEI]). 
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6 Appendix A 


Here we will show that, for the diverging flow, (i) there are no unstable modes if the basic flow 
is purely radial and (ii) all axisymmetric modes are stable. To do this, we employ the following 
theorem of Polya (problem 177 of Part V in [23], see also |24|). 

Polya’s theorem. Let the function f{t) he eontinuously differentiable and positive for 0 < t < 1, 
and also let fg f{t)dt exist. The entire function defined by the integral 

1 

J f{t)e^^dt = F{z) 

0 


has no zeros 


in the half-plane Rez > 0, if f'{t) > 0, 
in the half-plane Rez < 0, if f'{t) < 0. 


It should be noted that the interval (0,1) in the above theorem can be replaced by an arbitrary 
finite interval (a, 6 ). 

Consider first the case of purely radial flow. For 7 = 0 , the dispersion relation (1441) can be 
written as 

a 



<h(A:r) r dr 


(63) 


where 

^(s) = sInis)soK'^{sQ) - soI'^{sq)sKI^{s), So = ka. (64) 
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The change of variable of integration, ^ = r^/2, transforms (1631) to 

0^12 

D{a) = - [ (65) 

a J 
1/2 


If function f{^) were such that /(^) > 0 and f'{^) < 0 for ^ E then the above theorem 

implies that D{cr) has no zeros in the half-plane Rea > 0, i.e. all its zeros satisfy Rea < 0, which 
means that all modes are stable. 

The conditions for function /(^) that should be checked are equivalent to the following condi¬ 
tions for ‘h(s): 

<I>(s) > 0 and <0 for s E (A:, sq) • (66) 

It is convenient to introduce function T(s) by the formula 

T(s) = In{s)soK^{so) - solUso)Kn{s). (67) 


Conditions (j66p . expressed in terms of 'I'(s), become 

sT'(s) > 0 for sG{k,so), (68) 

(sT')^(s)<0 for sE(A:,so). (69) 


It is easy to see that function 'I'(s) is a solution of the modified Bessel differential equation 


1 A 

s ds 


ciT\ 


- - 


ds J 


1 + 


n 


T = 0 


and satisfies the following boundary conditions; 

T(so) = —1 and 'I''(so) = 0, 


(70) 


(71) 


where the first of these conditions follows from the Wronskian relation (1421) . 

First we prove the following auxiliary statement: 41 ( 5 ) 0 and 'I''(s) 0 for all A; < s < sq- 

To do this, we assume that either 'I'(s*) = 0 or T'(s*) = 0 for some s* E [A:,so)- Multiplying Eq. 
(I70l) by sTPs) and integrating from s* to sq, we find that 




/ 


s* 


T(s) {s^'Y {s)ds 


sT(s)^''(s) 




so 

S. 


Therefore, if either T(s*) = 0 or 'k'(s*) = 0, then 


J ^1 + ^‘^{s)sds = — J ^''^{s)sds, 


which is impossible. Therefore, both 'I'(s) and 'I''(s) must be nonzero for all A: < s < sq. 

Now we are ready to prove the required properties of 'I'(s). Since 'I'(so) < 0 and T(s) cannot 
change sign for s E [A:,so), we conclude that T(s) < 0 for all s E [A:, sq). Then, in view of the 
differential equation (fTOl) . we obtain 

(s^y (s) = s Tl -I- T(s) ^ (sT')'(s) < 0. 
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We have thus proved (j69|) . To prove ()68p . we observe that it follows from the differential equation 
(fTOl) and the boundary conditions that 'h"(so) < 0. This means that T'(s) > T'(so) = 0 at least 
near the end point s = sq- But since cannot change sign, it must be positive for all s G [k, sq), 
so that condition ([68|) is satisfied. 

Thus, the Polya theorem implies that for the purely radial converging basic flow (7 = 0 ), there 
are no unstable modes. 

It is easy to see that for the axisymmetric mode (n = 0) and for any 7, the dispersion relation 
(j44p also reduces to Eq. ()65p with n = 0, so that we may conclude that there are no growing 
axisymmetric modes. 


7 Appendix B 


Here we construct an asymptotic expansion of the solution to eigenvalue problem (I16p ~ (l26|) for 
large axial wave number k ^ 1. It is convenient to rewrite this problem in a form different from 
what has been obtained in Section 3.1.1. 

Let H{r) = ru{r). Then Eq. (|3ip can be written as 


in 

r 


n 


— H'{r) -[k^ + — ] G{r) = 0 


(72) 


where G = rv{r). It follows from (j72p that uJ 3 {r), given by Eq. (j22p . can be rewritten in term of 
H only as 


a 

UJ3 = — Or 

r 


Substituting this into Eq. (fTSl) and dropping the inessential factor in yields the equation 




in 


inj 0 ^ 

(J H - 2 - 1 - Or 


-dr 

r 


n- 


_j_ ^ 2^2 




Equation (fTSp must be solved subject to boundary conditions 
be written as 

i7(l) = 0, 77(a) = 0 

and either 

H'{1) = 0 

for the diverging flow (/? = 1) or 

H'{a) = 0 

for the converging flow (0 = —1). Equations (1751) and (|76l) follow from the incompressibility 
condition (l30p . 


= 0. (73) 

which, in terms of H, can 

(74) 

(75) 

(76) 


Diverging flow. Eigure [3] indicates that the azimuthal number of the most unstable mode behaves 
like n ^ k for large k. Therefore, in order to capture the stability boundary for large k, we consider 
the limit 

k 00, n —>■ 00, n ^ k. 

So, we set n = a fc in Eq. (I73P where a > 0 and does not depend on k. We also assume that 

7 = 7 A: and a = —i0a k^ + ak (77) 

where 7 = 0(1) and d' = 0(l)asA:—>-oo. Incorporating these assumptions into Eq. ([73]), we get 


1 A ~ 1 11^ 

I 3 - 1 ) + ^ + ^ - 


1 1 

r 


dr 


o? + 




= 0 . 


(78) 
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In the limit k ^ oo, this equation reduces to 


—i^a 




This implies that H (r) must be zero everywhere except a thin boundary layer near r = 1 where the 
above leading order term becomes small {0{k~^) as k ^ oo) and of the same order as some terms 
which we have discarded. To treat this boundary layer, we introduce the boundary layer variable 
^ such that 

r=l+i^ 

and rewrite Eq. ()78p in terms of At leading order, we obtain 


[a — 2ija ^ + 5^] 


1 




= 0 . 


(79) 


1 + 

Boundary conditions (1741) . (|75p . written in terms of take the form 

77(0) = 0, H'{0) = 0, 77 ^ 0 as I ^ oo. (80) 

The solution of Eq. (j79p , satisfying the first and the last of conditions (1801) , can be written as 


H{C) 


Cl 

2V2 


OO 

/ —as+i^as'^ 

0 




ds 


(81) 


where /i(a) = y/l + and Ci is an arbitrary constant. Note that formula () 8 ip is valid only for 
a satisfying the condition Re(cr) > 0. This means that our asymptotic result can only describe 
unstable eigenmodes. 

Substituting (fSTI) into the second boundary condition (|80p . we find that the condition of exis¬ 
tence of non-trivial solutions of problem (|79p . (1801) is 

OO 

7?2(cr,7,a) = J = 0 . (82) 

0 

Equation (f82|) represents the dispersion relation for eigenvalues a. Note that 7 ? 2 (d', 7 ,a), given 
by this formula, makes sense for a such that Re(cr) > and may have zeroes with —/i(a) < 

Re((T) < 0. However, only zeros of 7 l 2 (d', 7 , a) with Re((T) > 0 represent asymptotic approximations 
to the eigenvalues of the original problem. 

To make calculations easier, it is convenient to transform the dispersion relation to an equivalent 
form by deforming the path of integration on the complex plane of variables s from the positive 
real axis to the half-line: s = r G [0, 00 ). Then the dispersion relation takes the form 

OO 

0 

This equation can be further rewritten in term of the error function erf(z), but we do not do this 
here as Eq. (j83p is more convenient for numerical calculations. Typical roots of Eq. ()83l) are shown 
in Fig. [TOl Evidently, when 7 (a) is smaller than some critical value 7 cr(a), all roots are in the left 
half plane, and when 7 (a) > 7 ^( 0 ;), there is at least one root with Re(cr) > 0, which represents an 
asymptotic approximation of an unstable eigenvalue in the original problem. 
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Figure 10: Roots of Eq. (I83|) for a = 1 and 7 = 2, 2.5,4. 


To determine 7 cr(a) and Acr(a), we require a to be purely imaginary, i.e. a = iX with A E M. 
The dispersion relation (|83l) becomes 

OO 

0 

Equation (|84p has many roots. These roots corresponds to values of 7 (a) at which one of the zeros of 
the function D 2 {a, 7 , a) crosses the imaginary axis on the complex d-plane. We are only interested 
in the root that corresponds to the smallest value of 7 because it is this root that determines the 
stability boundary. In what follows we will consider only this root of Eq. (|84p . 

It turns out (the proof is below) that the function 7 cr(Q<) has a local minimum at a = 1 (i.e. 
when n = k as k ^ 00 ). Calculations yield 

min 7 cr-(a) = 7 cr(l) ~ 2.4671 and Acr(l) ~ 7.4331. 

a 

Since we are interested in the asymptotic behaviour of the stability boundary on the {k,'y), we 
choose a = 1 corresponding to the minimum value of Thus, the behaviour of the stability 

boundary in the limit /c —>• 00 is given by 


min 7 cr(a, k) = 2.4671 k + 0(1). 

This is shown by circles in Fig. [2j 

To prove that jcr(o:) attains its minimum value at a = 1, we change the variable of integration 
in Eq. (I84h : r = Qy/2/ix{a). This transforms (j84h to the equation 


where 


V2 

/i(a) 


CX) 

ye-7*O-eW4[*A*+7/2]C^^^0 

0 


1 * = 


2a 

/i2(a) 


7 , 


A* — 


V2 

ix{a) 


A. 


(85) 
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Then we make an observation that, up to an inessential constant factor, Eq. (|85p is exactly the 
same as Eq. (I 8 ^ for a = 1, with 7 and A replaced by 7 * and A*. This implies that, if 7 * and A* 
represent a root of Eq. ([85]), then 7 * = 7 cr(l) and A* = Acr(l)- This fact and the definition of 7 * 
and A* have a consequence that 


7cr(a) 


1 + 

2 a 


Tcr(l)) 


Acr(®) 


Vl + 


Acr(l). 


Since function /(a) = (1 + a^)/(2a) attains its minimum value at a = 1 and /(I) 
the required property of 7 cr(o)- 


1 , we obtain 


Converging flow. Eor the converging flow, a similar analysis shows that in the limit 


ka — 00 , n = a ka, 


the eigenvalue problem ([731) . (1711) . ([Tbil reduces to 


[a + 2ija 7 + dr,] 


1 

1 + a^ 




= 0 


( 86 ) 


and 


i7(0) = 0 , = 0 , 77 —>■ 0 as 7 —)• 00 . 

where rj, 7 and a are defined by 

T] = ka (^1 - ^ , 7 = 7 ka, ^ — [—ija{ka)'^ + a ka] . 


(87) 


It is easy to see that, the complex conjugate of Eq. ([86]) is equivalent to Eq. (17^ . with a replaced 
by its complex conjugate a. This means that if a and H{^) represent a solution of problem ([79p . 
([Ho]) , then a and and H{rj) solve problem (l86[) . ([871) . Therefore, the asymptotic result for the 
converging flow can be obtained from that for the diverging flow by simply replacing d by d and 
H{^) by Hence, we obtain 


min 7 cr(a) = 7 cr(l) ~ 2.4671, Acr(l) ~ —7.4331 

a 

where Xcr = Im(d) when Re(d) = 0. This means that in the limit ka —>■ 00 , n = a ka, we have 

min 7 cr(a, ka) = 2.4671 ka + 0(1). 


8 Appendix C 


Here we derive the asymptotic approximation ([58l) - ([62p . To obtain the regular part of the expansion 
(that is valid everywhere except the boundary layer near r = a), we substitute the asymptotic 
formula for the azimuthal velocity ([5T]) and Eqs. (l58]) - ([62]) into ([53[) ~ ([56]) . discard all boundary 
layer terms and collect term containing equal powers of 1 /i?. As a result, we obtain a sequence of 
equations, the first two of which can be written as 


KVq = 0, 

Rvj = -aivS + HvS. 

Here = (u^, for /c = 0,1 and operators K and B are defined as 


K^l = 


I + 

^ + ldr 


1 f.r _ 271 -r \ f) f,r \ 

-prVk+'^rPk 


+ + TPk 


y {^ + V^r)wl + ikpl 


( 88 ) 

(89) 


(90) 
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and 


/ Lul 
= Lvl 
\ Lwl 

In Eq. ([90]) . 

^2 

h{r) = ao — + inqi log r. 

and can be eliminated using the incompressibility condition 


dr (rug) + invg + ikrwg = 0. (92) 

Boundary conditions for Vq and are obtained by substituting (fMl) - ([6T]l into ([FT]) and collecting 
terms containing equal powers of 1/R. This yields 

^o(l) = ^o(a) = (93) 

i)S(l) = ^S(l) = 0, (94) 

i)S(a)+ i)o(0) = 0, Wo(a)+u;o(0) = 0, (95) 

■Ui(l)=0, {ti(a) + Uo(0) = 0, (96) 

DKl) = .hKl) = 0. (97) 


1 

:pi Uq 




1 -r I 2tn -r 
■piV0 + ^Ug 


( 91 ) 


We did not present boundary conditions for and wl at r = a, as they are not needed in what 
follows. 

Equation ([88l) and boundary conditions ([M|) and (fM|l represent the inviscid eigenvalue problem 
that was considered in Section 3. After it is solved, we know ag and Vq. Then boundary conditions 
(1951) are used to hnd the boundary layer corrections Vg and Wg. After that, the boundary layer 
part of the radial component of the velocity, Uq, can be found from the incompressibility condition. 
Once, Ug is known, Eqs. (IMl) and (f971) give us boundary conditions for Eq. (l89]l . which can then 
be solved, and the entire procedure can repeated as many times as necessary yielding higher order 
approximations. However, to find ui, we do not need to calculate the solution of (l8^ expicitly, all 
we need is to ensure that a solution does exist. Before describing how this can be done, we need 
to say a few words about the boundary layer. 

The boundary layer approximations are obtained as follows. We substitute Eqs. (f^ and (f58]l - 
([6^ into ([5^ - ([56]) and take into account that the regular part satisfies Eqs. ([88l) and ([89]l . Then 
we make the change of variable r = a(l — R~^ r]), expand every function of a(l — R~^ r]) in Taylor’s 
series at R~^ = 0 and, finally, collect terms of the equal powers of R~^. At leading order, the 
boundary layer equations are given by 


d‘^vl + dr^v^ = 0, 

(98) 

d^wl + dr^wl = 0, 

(99) 

—drjUg + invg + ikawg = 0. 

(100) 


The solutions of Eqs. (IM]l - (ll00h that satisfy boundary conditions (1^ and the condition of decay 
at inhnity are 


CX> 

Vg = —Vg{a) e“^, Wq = —Wg{a) e“^, Ug = — J (^invg{s) + ikawg{s)^ ds. 


V 
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Here the constant of integration has been chosen so as to guarantee that 'Uo(^) decays at infinity. 
Hence, the boundary condition for uKr) at r = a can be written as 

u[{a) = —{tQ(0) = —invQ{a) — ikawQ{a). (101) 

Now consider the non-homogeneous equation (|89l) . It has a solution only if its right hand side 
satisfies a certain solvability condition. To formulate it, we define the inner product 

a a 

(g, f) = y g-irdr = J {g^fi + 52/2 + ffa/a) rdr 
1 1 

where g = igi,g 2 , 93 ), f = ifi, f 2 , fs), and g is the complex conjugate of g. With respect to this 
inner product, we define the adjoint operator g 1 —?■ K*g by 

{g,K{) = {K*g,i) 


for any functions f and g satisfying the incompressibility conditions 


dr (rfi) + in f2 + ikrfs = 0 , dr (rgi) + in g2 + ikrgs = 0 
and the boundary conditions 

/i(l) = / 2 (l) = / 3 (l) = 0 , /i(a) = 0 , ( 102 ) 

51 ( 1 ) = 0, 51 (a) = 52 (a) = 53 ( 0 ) = 0. (103) 


Note that the boundary conditions for f and g are different. 

Now let g satisfy boundary conditions (I103p and be a solution of the equation 


^ - y 51 - ^ 5i + a ^ 


K*g = 


h (r) 


~ ^ dr] 92 + ^ §2 — ^ 5l + V “ 


^ - ^dr^ gs + ika 


= 0 


where function a can be eliminated using the incompressibility condition for g. Taking inner 
product of Eq. ([89]) with g gives us the required solvability condition: 


Qi — —o'iQ 2 + Qs 


where 

Hence, 


Qi = (g, Q 2 = (g, v);), Qs = (g, BvJ^) 

Qs — Qi 


fji = 


Q 2 


(104) 


Below are the explicit formulae for Qi, Q 2 and Q 3 that can be obtained after lengthy but standard 
calculations: 

a 

Qi = (n^ + k^a^) J e~^^^^Q{kr) rdr, 

1 

a 

<52 = —^ J e~^^^^^{kr) r^dr, 

1 

( a a 

J e~^^'"^W{r)^{kr) rdr + J {r"^ — Ij'^{kr) rdr ^ 
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where functions <h(s) and ^'(s) are defined in Appendix A and 


0(s) = Iniso)sK'^{s) - sl'^{s)Kniso) (sq = ka), 

W{r) = -fjo — + ( Y + (To(l - m7i) j (l + 7i) logr. 
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